Radar nonlinear multi-target tracking method with parallel PHD filter

Since probability hypothesis density (PHD) filters do not need explicit data association, they have recently been widely used in radar multi-target tracking (MTT). However, in existing PHD filters, sampling times are generally considered the same for all targets. Due to the limitation of antenna beam width in radar applications, the same sampling time for all targets will lead to a mismatch between the predicted data and measurement data, reducing the accuracy of radar MTT. In order to eliminate the estimation error with less computational cost, a radar nonlinear multi-target tracking method with a parallel PHD filter is proposed in this article. The measurement area is divided into several subspaces according to the beam width of the radar antenna, and the PHD of all subspaces is calculated in parallel. Then, multi-feature information in radar echo assists tracking and improves real-time performance. Experimental results in various scenarios illustrate that the proposed method can eliminate the estimation errors introduced by sampling time diversity at the cost of less computation cost, especially in cluttered environments.


Background The problem description of multi-target tracking in radar applications
In the practical application of radar, because of the limitation of antenna beam width, targets located in different spatial locations will be detected by different beams, resulting in sampling time diversity for different targets within the same scanning period.
To illustrate this problem further, Fig. 1 shows the surveillance area of the mechanical scanning radar; among them, the radar adopts the clockwise scanning mode, and the scanning speed is constant.In this scenario, there is a target moving at a constant speed.The target position is Z E k−1 at the end of the scan k − 1 .It will move from Z E k−1 to Z E k during scan k , and move from Z E k to Z E k+1 during scan k + 1 .Figure 1 also shows that the scanned positions of the target are respectively nearby Z 1 k , Z E k , and Z 1 k+1 .Assuming that the target has been stably tracked before the scan k , the target's predicted position by the filter will appear near Z 1 k in the scan k .Since Z 1 k is close to the target's predicted position, Z 1 k will be treated as the target-generated measurement, which is used to update the prediction priors.Similarly, measurement Z E k generated during scan k will be treated as clutter and does not contribute to the posterior calculation of the target.In the scan k + 1 , the target will move to Z E k+1 and be detected at Z 1 k+1 .According to the calculation result of the scan k , the target predicted position by the filter will be close to Z k+1|k .If the position is too different from the actual value, the filter may regard Z k+1|k as clutter and discard it.Meanwhile, the missed detection probability is used to compensate for the information loss due to the lost measurement, resulting in the divergence of the filter.Because the standard PHD filter does not consider the time information of the multi-target states, it is generally assumed that the measurements of targets are generated simultaneously, so the multi-target states cannot be accurately estimated.

RFS-based PHD filters
The PHD filter based on RFS is an approximate Bayesian filter that propagates the first moment rather than the entire multi-target posterior 10 .Among them, PHD is defined as an intensity function, expressed as 26 : The multi-target state at time k is naturally denoted as finite sets in RFS-based PHD filters: The multi-target measurement at time k is denoted by: where M(k) and N(k) are target number and measurement number respectively, F(X) and F(Z) respectively represent the collection of all finite subsets of the state space X and Z.
The PHD filter propagates the intensity function in time through two steps: prediction and update 10 .Let D k|k (x) denote the intensity function, the predictor operator D k|k−1 (x) is with b k|k−1 (x) designating the PHD of targets for spontaneous birth, b k|k−1 (x|x ′ ) designating the PHD of tar- gets for spawned, f k|k−1 (x|x ′ ) denoting the object transition kernel, and P S,k (x ′ ) denoting the target's survival probability.
The update operator is with g z (x) denoting the likelihood of individual targets, P D,k (x) denoting the target's detection probability, k c k (z) denoting the PHD of clutter, and

The parallel probability hypothesis density filter
As shown in Eqs. ( 4) and ( 6), the PHD filter does not consider the sampling time diversity of different targets, leading to the estimation error between prediction and measurement.Assuming that the radar needs N S beams to scan the whole measurement area, the measurement area is divided into N S independent subspaces.Since the scanning time of each subspace is very short, the sampling time diversity of each subspace can be ignored.Therefore, the multi-target state can be expressed as follows: with X S i denoting the multi-target state whose projection on the measurement space falls on subspace i , and X S i denoting the state subspace projected onto subspace i in the state space X .According to the above state model, the PHD can be expressed as follows: where symbol ⊎ stands for the disjoint union, for arbitrary sets A and B , C = A⊎B is equivalent to According to the independence of targets between subspaces, Eq. ( 9) can be further rewritten as the sum of PHD in different subspaces.An indicator function that accepts any arguments, including vectors, sets, etc., is denoted by www.nature.com/scientificreports/Then the PHD based on parallel processing is expressed as follows: with D S i (x) denoting the PHD of all targets in subspace i .Let D k−1 (x) denote the posterior PHD obtained at the k − 1 scan cycle.According to Eq. ( 11), the predicted PHD of the scan k is expressed as: with b S i k|k−1 (x) denoting the PHD of all birth targets in subspace i .The PHD of all birth targets in the surveil- lance area is The measurement model can be expressed as follows, which is analogous to the multi-target state model described in Eq. ( 8): where Z S i k is the measurement set from subspace i.The following is the written form of the PHD update operator: with g z (x) denoting the likelihood function of the measurement z .According to the independence of targets between subspaces, for arbitrary state x and measurement z , if x / ∈ X⊖ S i and z ∈ Z S i k , we have g z (x) = 0 .There- fore, Eq. ( 15) is further rewritten as: Vol.:(0123456789)The proposed method's specific steps In Eq. ( 16), PHD is denoted by the sum of the PHD components of different subspaces.Compared with the product form, this form of accumulation has better parallel processing ability, so it is suitable for parallel filtering scenarios in radar applications.This paper will use this framework to propose a radar nonlinear multi-target tracking method based on a parallel PHD filter that can eliminate the estimation errors introduced by sampling time diversity with less computation cost.The specific steps of the proposed method are shown in Fig. 2.

Measurement division
The measurement area is divided into N S subspaces; targets in the same subspace can be regarded as having the same sampling time.t S i k represents the sampling time of subspace i , the multi-target state is expressed as follows: with T S i denoting the sampling time space of the target in subspace i , and the multi-target measurement is expressed as follows: (16) The proposed method's specific steps.

Prediction
Each subspace takes the same sampling time, similar to the division of PHD in Eq. ( 11), then the PHD is expressed as: denote the posterior PHD of the scan k − 1 .According to Eq. ( 20), the predicted PHD of the k scan cycle can be defined by:

Feature-aided
The standard PHD filter uses Eq. ( 6) to update and iterate repeatedly.In the whole iteration process, g z (x) in Eq. ( 7) affects the filtering performance, but its calculation only uses the range and bearing information in the measurement.Considering that radar echoes contain abundant feature information in practical applications, such as range, bearing, amplitude, angle, Doppler frequency and SNR.Therefore, in this paper, Doppler and SNR contained in radar echo are considered as features to assist multi-target tracking because these two features of targets and clutter are different, and they are different for various targets.As a result, using these two features can reduce the relatively high computational cost by eliminating much clutter and distinguishing different targets.
In this paper, each target contains feature state Consequently, this paper uses feature-likelihood to refine the likelihood function of measurement in Eq. ( 15), as follows with g z (z k |x (i) k ) denoting the likelihood function of measurement z k generated by single target state x k of subspace i at time t S i k , p f (z (i) f k ) denoting the feature-likelihood of subspace i at time t S i k .Note that p f (z (i) ) and p f (z (i) P SNR ,k ) denote the Doppler frequency feature-likelihood and SNR feature-likelihood respectively, as follows: where σ 2 f d and σ 2 p SNR reflect the weight of the Doppler frequency feature-likelihood and SNR feature-likelihood in the calculation of the total likelihood function, in this paper, σ 2 f d = σ 2 p SNR = 0.85 .As a result, when measurements are closer to the clutter than to the target in probability, the likelihood function declines faster, and the clutter can be greatly reduced.In this way, the comparatively high computational cost is decreased while the real-time performance is enhanced.

Update
According to Eqs. ( 16) and ( 23), the PHD update operator at scan k can be defined by: with P D,t i k (x) denoting the detection probability, and t S i k c t S i k (z) denoting the clutter density at time t S i k .Finally, multi-target state estimation is extracted.This method divides the measurement area into several subspaces based on the radar antenna's beam width, and the PHD of all subspaces is calculated in parallel.At the same time, multi-feature information in radar echo is used to assist tracking.

Simulation and result analysis
Three simulation examples are used to assess the effectiveness of the proposed method.We compare the proposed parallel PHD filter based on feature-aided (FAP-PHD) against the TM-Joint-GLMB filter 24 , TM-PHD filter 25 , and the standard PHD filter through their Gaussian mixture implementation to simulate.The simulations in three different scenarios are shown as follows.

Simulation scene setting
Scenario 1: A mechanical scanning radar is located at (0 m, 0 m) and scans a semicircle region with the scenario setting range of [0 m, 2000 m] and the angle of [0°, 180°].The radar uses a bidirectionally continuous scanning mode, i.e., it scans repeatedly from 180° to 0° in a clockwise direction and then from 0° to 180° in an anticlockwise direction.The scanning speed is set to 180°/s.The surveillance region contains twelve targets.Figure 3 shows the true target trajectories with common scenes such as target birth, target death, target turning, and trajectory crossover, and Table 1 shows the initial states and lifetimes of targets.The semicircle region is equally divided into twenty subspaces; each subspace covers 9°.Radar range resolution is 10 m, and radar azimuth resolution is 1°; other scenarios are the same.In the surveillance region, clutter is evenly distributed, and the average clutter rate is = 50 points for each scan.In this scenario, the detection probability for each target is P D,k = 0.98, and the survival probability for each target is P S,k = 0.99.Scenario 2: Except for the clutter rate, all other parameters are identical to scenario 1, and we assess the efficiency of the proposed parallel PHD filter in this scenario with increasing clutter rates.Scenario 3: Except for the detection probability, this scenario is identical to scenario 1.We analyze the efficiency of the proposed parallel PHD filter in dense clutter ( = 300) with decreasing detection probability.

Target tracking setup
The single target motion model 27 adopted in the simulation is: www.nature.com/scientificreports/ In addition to the kinematic state, each target in the proposed FAP-PHD filter includes a feature state The nation ⊗ represents the Kronecker product, and the dimension of the identity matrix and a mean of 0, w k is Gaussian process noise.Here, the covariance matrix Q k|k−1 and the state-transition matrix F k|k−1 are respectively: with ϑ = 1 s denoting the maneuver correlation time, = 0.1 denoting the scalar acceleration standard deviation, and T s denoting the sampling period.
The model of single-target observation is where the observation matrix are the bearing and range of state x k , and e k is the observation noise.The observation noise e k follows Gaussian distribution with a covariance of R = diag([10, 1]) and a mean of 0. The Doppler frequency of each target in the proposed FAP-PHD filter is calculated by with f d and f r denoting the target's Doppler frequency and the radar frequency, c denoting the speed of light, and v r denoting the target radial velocity.The SNR of each target at adjacent moments follows a normal distribution with mean 6 and variance 2. For each scenario, the same target trajectory is used in 100 Monte Carlo trials to record the average performance, and each trial is independent of the other.Additionally, the efficiency of four filters is assessed using computational time and optimal subpattern assignment metric (OSPA) errors 28 .

Analysis of the results
Results analysis for scenario 1 Figure 4 presents the OSPA distance and its location and cardinality components for different filters.In this figure, "OSPA Dist, " "OSPA Loc, " and "OSPA Card" refer to OSPA distance, location components of OSPA, and cardinality components of OSPA, respectively.Figure 5 presents the computational time of different filters.It demonstrates that the FAP-PHD performs slightly better in location estimation and cardinality estimation than the TM-Joint-GLMB filter and outperforms the TM-PHD filter and the standard PHD filter, especially when targets disappear at k = 70 s and new targets are discovered at k = 20 s, k = 40 s, k = 60 s, and k = 80 s.However, the TM-Joint-GLMB filter is more computationally expensive, and the FAP-PHD filter performs marginally higher computation time than the standard PHD filter. (29) Table 1.The lifetimes and initial states of targets.When the clutter density is comparatively low, and the detection probability is relatively high, that is, = 50 and P D,k = 0.98, the OSPA error and computing time for different filters are shown in Table 2.When compared to the standard PHD filter, the proposed FAP-PHD filter reduces the average OSPA error by 66.88%, while the average OSPA errors of the TM-PHD and TM-Joint-GLMB filters are reduced by 19.23% and 51.42%, respectively.The average computational time of the FAP-PHD filter is increased by 1.4 times, while the average computational time of the TM-Joint-GLMB filter and the TM-PHD filter is increased by 18 times and 4.3 times, respectively, when compared to the standard PHD filter.As a result, the FAP-PHD filter enhances tracking performance compared to the TM-Joint-GLMB filter without significantly increasing time cost.

Results analysis for scenario 2
Figure 6 shows the mean OSPA performance improvement versus clutter rates.Table 3 shows the improvement in real-time performance versus clutter rates.As can be observed in Fig. 6, the average OSPA error of the proposed FAP-PHD filter remains essentially the same when the clutter density rises, but the mean OSPA errors of the other three filters increase.The time cost of the FAP-PHD filter is less than that of the other three filters, except that it is marginally higher than the standard PHD filter at = 50.
When the clutter density is relatively high, that is, = 300, the proposed FAP-PHD filter reduces the average OSPA error by 74.14% compared with the standard PHD filter.In comparison, the average OSPA errors of the TM-Joint-GLMB filter and the TM-PHD filter are reduced by 7.22% and 5.43%, respectively.The time cost of the FAP-PHD filter is reduced by 50.5%, while that of the TM-PHD and the TM-Joint-GLMB filters is increased by six times and 6.9 times, respectively.The results show that in this scenario, compared with the standard PHD filter, the tracking performance of the FAP-PHD filter is improved by 74.14%, and the real-time performance is enhanced by 50.5%.

Results analysis for scenario 3
Figure 7 shows the mean OSPA performance improvement versus detection probability.Table 4 shows the improvement in real-time performance versus detection probability.The results indicated that the average OSPA errors of all four filters increase with the decrease of detection probability.However, the proposed FAP-PHD filter obviously outperforms the TM-Joint-GLMB filter, the TM-PHD filter, and the standard PHD filter.Meanwhile, the average computational time of the FAP-PHD filter is lower than that of the other three filters.
When the clutter density is relatively high and the detection probability is relatively low, that is, = 300 and P D,k = 0.8, the average OSPA error of the FAP-PHD filter is reduced by 64.19% compared with the standard PHD filter.In comparison, the average OSPA errors of the TM-Joint-GLMB and the TM-PHD filters are reduced by  www.nature.com/scientificreports/5.82% and 3.81%, respectively.In the meantime, the time cost of the FAP-PHD filter is reduced by 55.1%, while that of the TM-PHD and the TM-Joint-GLMB filters is increased by 5.3 times and 5.8 times, respectively.The results show that in this scenario, compared with the standard PHD filter, the FAP-PHD filter has a tracking performance improvement of 64.19% and a real-time performance improvement of 55.1%.All simulation results in three scenarios demonstrated that the proposed FAP-PHD filter performs exceptionally well in real-time and tracking, especially in cluttered environments.The proposed method adopts the parallel probability hypothesis density filter, which divides the measurement area into several subspaces and then takes the SNR and Doppler frequency contained in radar echoes as the features of auxiliary multi-target tracking.These features can eliminate clutter and differentiate targets from each other, thus eliminating the estimation errors introduced by sampling time diversity with less computation cost.As expected, with increasing clutter rate and decreasing detection probability, the FAP-PHD filter has a lower mean OSPA error and a faster computing time when compared to the TM-Joint-GLMB filter, the TM-PHD filter, and the standard PHD filter.Simulations are consistent with the theoretical analysis.
In summary, simulations demonstrated that the proposed method not only outperforms the TM-Joint-GLMB filter, the TM-PHD filter, and the standard PHD filter in tracking performance but also has advantages in realtime performance, especially in cluttered environments.Therefore, the FAP-PHD filter is suitable for radar application scenarios that require high real-time and tracking accuracy performance.

Conclusion
This article proposes a radar nonlinear MTT method with a parallel PHD filter.According to the beam width of the radar antenna, the measurement area is divided into several subspaces to realize the parallel filtering on each subspace.Then, the feature information in the radar echo is used to eliminate clutter further, thus eliminating the estimation errors caused by sampling time diversity with less computation cost.Simulations demonstrated that the proposed method could address the problem of sampling time diversity in radar applications and obtain higher tracking performance with less computational cost, especially in cluttered environments.However, the proposed method still has challenges, such as extended target tracking and multi-sensor fusion.In further research, we focus on applying the proposed method to multi-sensor fusion and extended target tracking.

Figure 1 .
Figure 1.Surveillance area of the mechanical scanning radar.
where x f d ,k represents the Doppler frequency, x P SNR ,k represents the SNR of the target at time k , [ẋ k , ẏk ] T and [x k , y k ] T represent the velocity and the position, respectively, [ẍ k , ÿk ] T represents the acceleration.The state of each target is the concatenation of kinematic and feature state, denoted as xk = [x k , x f k ] T .Similarly, the corresponding target observation zk = [z k , z f k ] T is the concatenation of kinematic and feature measurement, where kinematic measurement

Figure 6 .
Figure 6.The mean OSPA performance improvement versus clutter rates.

Figure 7 .
Figure 7.The mean OSPA performance improvement versus detection probability.

Filter P D,
k = 0.8 P D,k = 0.85 P D,k = 0.9 P D,k = 0.95 P D,k = 0.98 δX S i

Table 2 .
OSPA error and computing time for different filters ( = 50, P D,k = 0.98).Significant values are in bold.

Table 3 .
The improvement in real-time performance versus clutter rates ( P D,k = 0.98).Significant values are in bold.

Table 4 .
The improvement in real-time performance versus detection probability.